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We present a first principles theoretical analysis of the entanglement of two superconducting 
qubits in spatially separated microwave cavities by a sequential (cascaded) probe of the two cavities 
with a coherent mode, that provides a full characterization of both the continuous measurement 
induced dynamics and the entanglement generation. We use the SLH formalism to derive the full 
quantum master equation for the coupled qubits and cavities system, within the rotating wave and 
dispersive approximations, and conditioned equations for the cavity fields. We then develop effective 
stochastic master equations for the dynamics of the qubit system in both a polaronic reference 
frame and a reduced representation within the laboratory frame. We compare simulations with and 
analyze tradeoffs between these two representations, including the onset of a non-Markovian regime 
for simulations in the reduced representation. We provide conditions for ensuring persistence of 
entanglement and show that using shaped pulses enables these conditions to be met at all times 
under general experimental conditions. The resulting entanglement is shown to be robust with 
respect to measurement imperfections and loss channels. We also study the effects of qubit driving 
and relaxation dynamics during a weak measurement, as a prelude to modeling measurement-based 
feedback control in this cascaded system. 


I. INTRODUCTION 

Entanglement between remote parties is a key resource 
in many quantum information applications, including 
quantum teleportation, quantum key distribution, and 
quantum metrology, and is also central to the notion 
of a quantum network [1]. Distributed entangled states 
are also a critical component for scalable quantum com¬ 
puting, since they enable long-range gates between spa¬ 
tially separated qubits [2]. Accordingly many different 
approaches have been proposed to distribute or gener¬ 
ate entangled states among systems that are significantly 
spatially separated. Distributing entangled states af¬ 
ter preparation at a central location is practically chal¬ 
lenging since decoherence in distribution channels typi¬ 
cally degrades entanglement, e.g., [3, 4]. Alternatively, a 
long-range coupling between remote systems can be en¬ 
gineered by exchanging single quanta, and entanglement 
can be generated this way, as has been recently demon¬ 
strated for atoms, photons, and combinations thereof 
[5, G]. 

A fundamentally distinct approach for preparing en¬ 
tangled states of systems residing at remote locations is 
to perform a, joint measurement on them. Most proposals 
for achieving such joint-measurement-enabled entangle¬ 
ment interfere photons that are spontaneously emitted 
by atoms (or artificial atoms) in such a way that sub¬ 
sequent detection of a photon makes the identity of the 
emitter indiscernible, e.g., [7-12], and thus projects the 
remote atoms into an entangled state. The degree of en¬ 
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tanglement generated is heavily dependent on both the 
quality and stability of the interferometer and efficiency 
of detection of spontaneously emitted photons. As a re¬ 
sult achieving high fidelity entangled states with this ap¬ 
proach is challenging, although several proof-of-principle 
experiments have demonstrated validity of the approach 
[13, 14]. An alternate approach is to perform a joint mea¬ 
surement by sequentially interacting two systems with a 
coherent light mode. This has been explored as a method 
for generating entanglement theoretically [15, 16] and ex¬ 
perimentally implemented using collective excitations of 
atomic clouds [17]. Most recently, a sequential probe 
has been utilized to probabilistically entangle supercon¬ 
ducting qubits in separate microwave cavities [18]. In 
this work we develop a theoretical description of that ex¬ 
periment from first-principles, providing a rigorous and 
general theoretical framework for the generation of en¬ 
tanglement by joint dispersive measurement of qubits in 
distinct cavities and analyzing in detail the potential and 
limitations of entanglement generation in this setting. 

In the dispersive interaction regime, a coherent mode 
reflected off a cavity with an embedded qubit acquires 
a phase shift that depends on the internal state of the 
qubit. This motivates the essential idea behind the en¬ 
tanglement generation scheme we study here, namely, to 
perform a measurement of the parity of the qubit pair 
excitation state by sequentially probing the two cavities 
that contain them, and performing a homodyne measure¬ 
ment of the total phase acquired by the twice-reflected 
probe field. Fig. 1 shows a schematic of the appara¬ 
tus. Ideally, the qubit observable that corresponds to 
this measurement, which we shall refer to as a half¬ 
parity measurement, takes the form Otp = cr\ + cr\, 
where cr\ is the Pauli-Z operator on the qubit. This 
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observable cannot distinguish between the qubit basis 
states 101) and |10) Therefore if the initial state of 
the two qubits is the equal superposition state |'I'o) = 
|(| 00 ) + | 01 ) + | 10 ) + | 11 )), the ideal half-parity mea¬ 
surement will yield the states | 00 ) or | 11 ), each with a 
probability 1/4, or the state :^(|01) -I- |10)) with proba¬ 
bility 1/2. This is to be distinguished from the full par¬ 
ity measurement oi Ofp = cr),cr^, which yields the states 
;^(|00) -I- |11)) and ^(|01) -I- |10)) with equal proba¬ 
bilities. In the following we develop a detailed model of 
this sequential probe measurement from first principles, 
including all non-idealities present in the experiment of 
Ref. [18]. Although we develop the model within the 
context of the superconducting qubit experiment of Ref. 
[18], it applies more generally to any implementation of 
cavity-QED, including in the optical domain. 

The remainder of the paper is structured as follows. 
Section IIA presents the dynamical model for full sys¬ 
tem derived from cascaded systems theory [19-21]. Then 
in section IIB we perform a two-cavity polaron trans¬ 
formation on the model to obtain an exact, dressed de¬ 
scription for the qubits alone that is easier to simulate 
than the full dynamical model. This polaron transform 
and subsequent derivation of an effective master equa¬ 
tion for the dressed qubits constitutes a generalization of 
the methods first presented in Ref. [22] . In Sections IIC 
and IID we simulate the qubit-reduced dynamics in the 
laboratory frame, including the loss and revival of coher¬ 
ence between qubits. Section III derives physical require¬ 
ments and criteria for generating entanglement between 
the remote qubits. Section IV provides simulation data 
for a range of realistic experimental parameters and dis¬ 
cusses the viability of obtaining high-grade concurrence 
betwen the qubits. Section V develops a perturbative 
treatment of qubit driving during the continuous mea¬ 
surement. Finally, section VI provides a summary and 
assessment of the benefits and possible extensions of this 
approach for other quantum processing tasks with super¬ 
conducting qubits. 


II. DERIVATION OF THEORETICAL MODEL 

A. Full model of cascaded cavities 

Consider the apparatus shown in Fig. 1. Fach cavity 
has two ports, with asymmetric transmittivities. The 
“input” port on each cavity is low transmittivity ( 7 ^) and 
the “output” port is high transmittivity (k^). The probe 
field e{t) interfaces with the output ports of both cavities 
a distance L apart (in [18] a distance L = 1.3m was 
achieved), before impinging on the homodyne detector. 


^ Here, and in the following, for conciseness we omit tensor prod¬ 
ucts when writing multiparty states. 



FIG. 1: Sequential probe of two spatially separated cavities, 
each containing a qubit coupled dispersively with strength Xi 
to its fundamental cavity mode. The cavities are asymmetric 
with the K ports being more transmissive than the 7 ports. 
The beam-splitter between the cavities models the losses in¬ 
duced by the circulator that enforces one-way field propaga¬ 
tion between cavities. We allow for arbitrary coherent state 
drives {Ad{t) and Bd{t)) into the weakly coupled ports of both 
cavities. The output field that results after sequential reflec¬ 
tion of the probe field from both cavities, z{t), is measured 
by a homodyne detector with efficiency 7 m and at a phase (j) 
with respect to the probe e(t). 

The cavities are operated in the dispersive regime, where 
the Hamiltonians in the two cavities are given by 

Ha = Aia'^a-I-xia'^acr/, 

Hb = A2b^b + X2b^bcr^, (1) 

respectively, where a(b) is the annihilation operator for 
the fundamental mode in cavity 1(2), is the Pauli 0 
operator for qubit 1 ( 2 ), Ai = ojd — w/ is the detuning of 
cavity i from the probe field e(t) (frequency ujd), and Xi is 
the qubit-cavity coupling in the dispersive regime. These 
Hamiltonians are in the interaction frame with respect 
to the free Hamiltonians for the qubits: —^crl — 

The input ports for the two cavities can also be used 
for driving the cavity or qubits at their respective fre¬ 
quencies (with fields Adit) and Bd{t), respectively), for 
state initialization and tomography. We will see that 
the coherent drive Bd{t) will be useful for compensating 
against asymmetries in the parameters between the two 
cavities and qubits. To minimize back-reflection of the 
probe held and ensure its unidirectionality, a circulator 
is inserted between the two cavities. Losses associated to 
this circulator will be included in the model developed 
below. The parameters of the system, including cavity 
transmitivities, losses, cavity coupling, are labeled in Fig. 
1 . 

The dynamical model for the apparatus described in 
Fig. 1 can be derived using the cascaded cavity theory 
of Gardiner and Carmichael [19, 20], or by the modern 
SLH quantum network theory [21, 23]. We utilize the 
latter here and Fig. 2 presents an SLH network diagram 
that is equivalent to the apparatus in Fig. 1. Fach block 
Gi is specified by an {S,L,H) triple, where S' is a scat¬ 
tering matrix, L a coupling matrix and H a self-energy 
matrix. Gi , G 2 and G 3 represent coherent displacements 
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of the input vacua, G 4 and Gg represent the cavity-qubit 
systems, and G 5 represents a beamsplitter modeling the 
lossy circulator. The output held z{t) emerges from the 
output port of cavity 2 and is monitored by homodyne de¬ 
tection (see below). The SLH triples {S, L, H) [21, 23, 24] 
for these blocks are 


Gi 

G 2 

G 3 


Ga 


Gs 

Gg 


(l,e-(t), 0 ) 


-1 


, Aia'l'i 


’ L . 

j/m Vi - vi 
- m y/m 


-f Xia^ao-i 

, 0 , 0 ) 


— I*; 




, A 2 b^b + X2b^bcr) 



(a) 



where I 2 is a 2 x 2 identity matrix and rn is the efficiency 
of the circulator between the cavities {i.e., rji = 1 implies 
no loss) The SLH representation of the entire system 
is then formed by performing the following concatenation 
(ffl) and series (<l) products [23]: 

(Go ffl G^^^ ffl Go ffl G^^^) < (Gs ffl Go ffl Go) < 

(GofflG4fflGo)<(GofflGifflG2fflG3), (2) 

where Go = ( 1 , 0 , 0 ) is a pass-through component, and 
we have split the two ports of the second cavity as 


(b) 

FIG. 2: SLH network decomposition of the apparatus in Fig. 
1. The SLH triples {S,L,H) for each block are specified in 
the main text, (a) shows an SLH block representation super¬ 
imposed on the experimental apparatus, and (b) shows the 
SLH block diagram redrawn more conventionally with inputs 
on the left and outputs on the right. The explicit forms for 
the resulting S (scattering), H (self-energy), and L (coupling) 
matrices for the overall model are given in the Appendix. All 
inputs are in the vacuum state (indicated by dotted lines in 
(a)) and the single monitored output is z{t), the field reflected 
from the output port of cavity 2. All other outputs are not 
monitored: this is indicated in (a) and (b) by their termina¬ 
tion. 


Gg^' = (^-1, v^b, A2b^b-hX2b^bcr2^ 

= (-l,v^b, 0 ) 

for convenience (without this splitting, we would have 
to insert a routing element to swap the third and fourth 
signal lines after G 4 in Fig. 2). The key assumption in 
this SLH representation of the entire network in terms of 
its components is that the fields propagate with negligible 
time delay between the components, which we assume to 
be true. 

Evaluating the series and concatenation products in 
Eq. (2) using the rules specified in Refs. [21, 23] yields 
the overall G = {S,L,H) for the network, from which 
an equation of motion for the two qubits and inter-cavity 
modes may be extracted (see Appendix). Adding phe¬ 
nomenological Markovian dephasing terms for the qubits 
with rate parameters 7 ^,* = 1 , 2 , then results in the fol¬ 
lowing master equation for the cavity mode and qubit 


^ The zero (0) elements in SLH triples should be interpreted as 
zero matrices or vectors of the appropriate dimension. 


degrees of freedom: 

p] -I- CcQ + CqQ 

'D[yjKi{l - rii)a] + 7iX>[a]p -h 72T>[b]p 
-l-T>[-^Ki?7;a -I- v^b]p 
2 

(3) 

i=l 

where g is the combined density matrix of the two cav¬ 
ity modes and qubits, and I?[A]H = ABA^ — AB — 
\bA^A. The effective Hamiltonian for the coupled sys¬ 
tem is 

H' = Ha + Hb + K + K 

K = -f'^(a^b-b^a), 

H'^ = t(A,(t)at-A5(t)a + i?,(t)bt-H2(i)b),(4) 

where K 12 = and the effective cavity drives are 

Ad(t) = y/^Ad(t) + 

Bd{t) = \/^Bd(t) - (5) 


dt 

CcQ = 

BqQ = 
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iJ' describes the effective direct coupling of the two cav¬ 
ity modes due to the probe field that interacts with both 
cavities. Similarly, the effective drives for the cavity 
modes in are composed of the probe field e{t) and 
the drive fields entering the input ports. The change in 
the sign of e{t) in the two expressions for effective cavity 
drives reflects the phase shift that the probe field picks 
up as it reflects from cavity 1. Note that the coupling 
in is reduced by the factor due to losses in the 
circulator. Although the Hamiltonian form of this cav¬ 
ity coupling looks reversible, the irreversibility enforced 
by the circulator is nevertheless captured in the network 
model when the dissipative dynamics modeled by Cc is 
included. We shall see the effects of this explicitly below 
(see also Ref. [20]). 

Cq accounts for intrinsic decoherence for the qubits, 
which is assumed to be pure dephasing dynamics. In 
most of this work we neglect the contribution of re¬ 
laxation (and heating) which typically contributes on a 
timescale Ti much longer then the timescales of interest 
for establishing entanglement by continuous joint mea¬ 
surement, e.g., [18]. However, in section V we discuss 
how this effect may be included in the reduced model for 
the coupled qubit dynamics in the laboratory frame that 
is derived in section HD. 

Each dissipator term in Cc takes into account the ef¬ 
fect of a field irreversibly coupling out of the combined 
system. In particular, T>\\J ki(1 — rii)a] accounts for pho¬ 
tons lost between the two cavities, I?[a] accounts for the 
field emitted from the input port of cavity 1, I?[b] ac¬ 
counts for the field emitted from the input port of cavity 
2, and T>[—ydtirjla -|- yd«2b] accounts for the probe field 
that is sequentially reflected off the output ports of cav¬ 
ities 1 and 2. This final output channel is the only one 
that is monitored, and this correlated dissipator encodes 
the fact that when coherent light escapes the system, 
it has interacted with both the first and second cavity. 
It is therefore impossible to distinguish which cavity a 
decayed photon has come from, and thus it must be de¬ 
scribed using a combined operator. To complete the full 
model in the presence of measurement, we must describe 
the evolution of the system conditioned on homodyne 
measurement of the output field from cavity 2, 

z{t) = + yK2b(t). (6) 

The homodyne measurement is implemented by mixing 
this signal field with a local oscillator of fixed phase ref¬ 
erence, (j), with respect to the initial phase of the probe 
field. This phase reference sets the measurement quadra¬ 
ture. The corresponding time evolution is given by [25] 

— = —i[H\g\+Ccg + Cqg-\-CmQ 

Cmg = + ^/K^b)]g, (7) 

where 0 < ?7m < I is the efficiency of the measurement, 
(j) defines the measurement quadrature, and ^(t) is Gaus¬ 
sian white noise due to the measurement. This equation 


is in Ito form [26] and therefore = dW{t), where 

dW{t) is a Wiener increment satisfying E{dW{t)} = 0 
and E{dW(t)dW{s)} = d{t — s) {E denotes expectation 
value). The nonlinear conditioning superoperator "H is 
defined as: n[A]B = AB + BA^ - Tt{AB + BA^)B. 
Eq. (7) describes the conditioned state of the system un¬ 
der a homodyne measurement trace of the voltage 

V{t) = + V^b)) + ^(t), (8) 

where (A) = tr(Agi). Eq. (8) expresses the monitored 
voltage in terms of the measured observable z{t), which 
is a linear combination of intra-cavity field operators a 
and b (Eq. (6)). 

The Heisenberg equations of motion for the expected 
values of the intra-cavity fields under the unconditioned 
evolution described by Eqs. (3)-(4) are: 

(a) = -tAi(a) -ixi(o-^a) - {a) + Ad{t) 

(b) = -»A2(b) - ix2{a-lb) + Ki2(a) _ (b) 

+ Bd{t). 

These evolution equations make explicit the fact that 
the second cavity (b) is driven by the first, (a), but 
not vice-versa (i.e., the irreversibility of the coupling 
between cavities). We assume that the driving fields 
{e{t), Ad{t), Bd{t)) are all coherent states and therefore 
these expectation values are simply the coherent state 
amplitudes of the intra-cavity fields. We can write these 
coherent state amplitudes conditioned on the qubits be¬ 
ing in specific states as: 

= -zAiAM _ + Tl ^(r) 

+ Ad{t) 

+Bd{t) + ki2 A^'^\ (9) 

Here A = {a),B — (b) and the superscripts r, s G {0,1} 
indicate the conditioning on the state of the first and 
second qubit, respectively. The state of the second cav¬ 
ity, B^^^\ is conditioned on the states of both qubits but 
the state of the first cavity, A^A ^ jg only conditioned on 
the state of the first qubit, since there is no information 
flowing back from the second to the the first cavity. Ex¬ 
plicitly, = AA) and 

These conditioned equations for intra-cavity amplitudes 
are linear and can be solved exactly (first solving for A 
and then for B) for any values of the driving fields. Their 
exact solutions will be used below. 


B. Dynamics in the polaron frame 

While the conditioned dynamical equation in Eq. (7) 
is a full model of the experimental setup in Ref. [18], it 
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is difficult to simulate since it involves both qubit and 
cavity degrees of freedom. Therefore it is convenient to 
derive an effective SME for the qubit degrees of freedom 
only. For this purpose, in this section, we develop an 
SME in a polaron frame where the average state of both 
intra-cavity fields is displaced to the vacuum. This SME 
is exact within the rotating wave approximation (RWA) 
and dispersive approximation implicit in Eq. (1), and be¬ 
comes easy to simulate since the intra-cavity fields are 
always in the vacuum state. 

The polaron transformation provides a representation 
in which the cavity and qubit degrees of freedom are hy¬ 
bridized. The correct transform in this two cavity case is 
g^{t) = U{ty Q{t)U(t), with 

m = E.., , (10) 

where V\ij = (i| ® (j| projectors onto qubit 

states and D^^yX] is a displacement operator for cavity 
field 1(2), i.e.. 


Performing all the polaron transformations, the SME de¬ 
scribing conditioned evolution in this frame then takes 
the explicit form (for an unnormalized density matrix): 

do^ 

= —i\HA + Hb + H'^ Hq, g^] 

+Ccg^ + C,qg^ + Cjqg^ + ^mg^ + C,mqg^ 
+a[g^, TiflJj -|- «:i2n|,] -|- [riFlo -f Ki2r\b, 
r2n]^ -1- «;i2n)j] -t- [r2nb -y «;i2na, 

(14) 

where 

Hq = t(^Ayt)ni-A*^{t)na + Byt)nl-B:{t)nb) 
^'q = 'y 2 'B[r\b] + {ki{ 1 - r]l) + "f^Vlna] 
+'D[-^Kir]ir\a + v^rih] 

Bmq = y/V^C{t)n[e^‘^{-y^KAmr\a + (15) 


Di[X] = 

D2[X] = 

For convenience, we define the following time-dependent 
qubit operators that depend on the intra-cavity field 
states: 

n,(t) = noA(o)(t)-f niA(i)(t), 

nfe(t) = ^ n,,R(*^)(f), (11) 

i,1=0,1 

where Flo = floo + rioi and Hi = flio-l-nii. The quanti¬ 
ties AW(t),R(’'*)(t),n a/b{t) are all time-dependent, how¬ 
ever, in the following we will often suppress the time pa¬ 
rameter for the sake of brevity. Eq. (11) can also be 
viewed as qubit projectors whose relative amplitudes de¬ 
pend on the cavity fields, which arises as a consequence 
of the hybridization of cavity and qubit degrees of free¬ 
dom induced by the polaron transform. Note that with 
this definition, the polaron transformation can also be 
written as 


U{t) = Dq[na{t)]D2[nb{t)] (12) 

The temporal evolution of the joint density matrix in 
this polaron frame is given by 

= -^[H'P,gP]+CygP + CygP + CPgP 
-U’^Ug^ - g^il'^U, 

where H'^, are as defined in Eq. (3) and 

Eq. (7), but with each operator transformed into the po¬ 
laron frame according to = U{tyOU{t). For exam¬ 
ple, transforming the field annihilation operators yields 

U{tyaU{t) = a + no(t) 

U{tybU{t) = b + nb(t). (13) 


with Fi = 'ji + Ki being the total decay rate of cavity 
i, and 'H[A]i3 = AB + BA^. Comparison with Eq. (7) 
shows that the polaron transformation has resulted in the 
additional terms Hq,Cq and C^q in the master equation. 
We sacrifice normalization in the following for simplicity, 
noting that the normalizing factor can always be recov¬ 
ered by computing tr (qf^)- In this polaron frame, there 
is no drive of the cavity modes by Ad and Bd, because 
we are dynamically shifting the cavity states back to the 
vacuum. As a result, if the cavities are unpopulated ini¬ 
tially, they remain unpopulated at all times. One way to 
see this is to note that all field operators act as annihila¬ 
tion operators on g^ in Eq. (14), and therefore there is no 
change in the states of the two cavity modes if they start 
in the vacuum. As a consequence, all terms involving 
cavity mode operators a, b in Eq. (14) have no effect and 
can be dropped. This results in the following equation of 
motion for the qubit degrees of freedom in the polaron 
frame: 


dg^ 

dt 


-i[Hq, g^] + Cqg^ + Cqg^ + Cmqg^ (16) 


The terms in Hq represent Stark shifting of the energy 
levels due to the interaction with the measurement pulse. 
The terms in Cq represent information about the qubits 
leaking out (and thereby also dephasing the qubits) as 
a result of light exiting the various output ports in the 
system. Cmq represents stochastic measurement noise on 
the system as a result of the monitoring of one of the 
output ports. 

It is important to note that Eq. (16) contains no ad¬ 
ditional approximations beyond the RWA and the ap¬ 
proximation of the Jaynes-Cummings interaction by the 
dispersive interaction between qubits and cavity modes. 
However, due to the polaron frame transformation, as 
long as the cavities are initially unpopulated, this equa¬ 
tion efficiently simulates the coupled qubit and cavity 
degrees of freedom without the cost of keeping track of 
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the quantized field states in the cavity. Instead, their in¬ 
fluence is captured by the time-dependent operators 
and rift in Eq. (16). 


C. Transforming back to the lab frame 

In order to make predictions with respect to the lab 
frame, we transform the density matrix that results from 
Eq. (16) back into the lab frame. We can achieve this by 
first noting that the state of the system at an arbitrary 
time in the polaron frame takes the form: 

Q^{t) = \ij) {kl\ (g) |00) (00| , (17) 

ijkl 

where rijki{t) is the solution to Eq. (16). The ijkl in¬ 
dices run over {0,1} and index the qubit states. The 
second term in the tensor product, |00) (00|, is the state 
of the intra-cavity fields. Both modes are in the vacuum 
state since in the polaron frame the cavities remain un¬ 
occupied. Then the state of the entire system in the lab 
frame is given by Q{t) = U(t){t)W (t), and the state of 
the qubits in the lab frame is given by 

p{t) = trci,c2 {U{t)g^{t)U\t)) , (18) 


D. Reduced equation of motion for the qubits 

Another approach for obtaining the state of the two 
qubits at any time is to formulate an equation of motion 
for just the qubit degrees of freedom in the lab frame. 
We begin with the expression for a general two qubit 
state in the lab frame given in Eq. (19). Taking the time 
derivative of this state yields: 

for the diagonal elements, and 

-^Pijkl = ® ^ + Pijkl 

( 22 ) 

The time derivative of rijj-i is determined by the polaron 
frame SME of Eq. (16), and the time derivative of the 
compensation factor can easily be found from its defini¬ 
tion in Eq. (21) together with the equation of motion for 
the conditional intra-cavity fields, Eq. (9). This is simi¬ 
lar to the approach taken in Ref. [22] for a single cavity 
setup. 

Computing the derivatives required in Eq. (22) and 
canceling common factors yields the following equations 
of motion for the components of the (unnormalized) lab 
frame qubit density matrix: 


where the trace is taken over both cavity modes. Writing 
P{t) = Pijkiit) \ij) {kl \, (19) 

ijkl 


we can determine the relation between pijki{t) and 
’I’ijkii't)-, using the definition of the polaron transform in 
Eq. (10) and evaluating Eq. (18). We find that the diag¬ 
onal elements remain unchanged by the transformation, 

Pijijit) = 

but that the off-diagonal components are modified as 

Pijkiit) = (20) 


where the compensation factor is 


P^Jkl = p^Jkl[^2xl{l-S,k)[i-lyA^'^'>*A(■^'>^ 


with 


-27d(l-<^*fc) 

-2xj{l-Sji) 



(23) 



(24) 


being the conditional output fields. This evolution can 
also be represented in matrix form (again for an unnor¬ 
malized qubit density matrix) as: 

^ = '^aijkiit)r\^jp{t)r\ki + Cqp{t) + CjnqP{tX25) 

ijkl 


1^4(0 _ yl(fc)|2 |5(b)_^(fe0| 


■ ( 21 ) 


These relations suggest that an efficient method for 
simulation of the system in the lab frame is to com¬ 
pute the time dynamics in the polaron frame according 
to Eq. (16), giving rijki{t), and then to compute the com¬ 
pensation to the off-diagonal elements given by Eq. (20) 
at each time, to get the reduced state of the qubits in the 
lab frame. 


with 


+12x2(1-Sji) 

Eq. (25) should be treated with care since although it 
looks like an SME with the deterministic component in 
Lindblad form, it is not strictly in Lindblad form because 
the coefficient matrix defined by aijki{t) is not neces¬ 
sarily positive. Hence this equation can result in non- 
Markovian evolution of the qubit. Physically, this is a 
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FIG. 3: (Color online) Distribution of qubit state populations 
(200,000 trajectory simulations) as a function of normalized 
homodyne voltage and measurement pulse width, for simula¬ 
tion parameters corresponding to slightly different qubits. See 
section IIE for the parameter values. The initial state is given 
in Eq. (26). The z-axis is the total number of trajectories that 
result in the state |00) (brown), |01)(blue), |10)(green), or |11) 
(yellow). 
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FIG. 4: (Color online) Bottom panel: Sample trajectory of 
|poiio(t)|, the absolute value of an off-diagonal qubit density 
matrix element simulated using the three dynamical equa¬ 
tions derived in section IIB (polaron, Eq. (16)), IIC (lab 
frame compensated polaron, Eqns. (18)-(21), and IID (re¬ 
duced equation of motion, Eq. (25)). We use the ideal setting 
with identical cavities and lossless transmission. All parame¬ 
ters are described in the main text. Top panel: measurement 
pulse Adit) for this trajectory. 


result of tracing out the intracavity degrees of freedom 
that are strongly entangled with the qubit states. There¬ 
fore physical interpretation of the rate coefficients aijkiit) 
in Eq. (25) is difficult. Nevertheless, Eq. (25) does gen¬ 
erates the correct qubit evolution in the lab frame and 
presents an alternative to simulating the qubit evolution 
in the polaron frame according to Eq. (16). 


E. Simulation in different frames 

Simulation of the system dynamics using any one of 
these three dynamical equations allows the generation of 
statistics for the homodyne voltage measured in individ¬ 
ual runs of an experiment. As discussed above, in both 
the polaron and laboratory frames (Eqs. (16) and (25)), 
the measurement observable corresponds to 

Re{e'’‘^{-y/Kir]ina + v^Hf,)}. 

Fig. 3 demonstrates the distribution of measurement 
outcomes as a function of the normalized homodyne 
voltage and the measurement pulse width. We used 
the polaron frame reduced master equation, Eq. (16), 
with a large number of random realizations of dW (t) 
(200,000 trajectories) to obtain these results. This choice 
of frame is made here for computational efficiency. For 
these simulations we use the parameters ^ = 2MHz, 
g = l.OMHz, = 18MHz, = 22MHz, a" = A 2 = 0, 
71 = ioA2 = io’ = 13MHz, 77/ = 0.9, and 

bm = 0-4. Given an initial state that is the equal super¬ 
position of the computational basis states, 

|4^o) = ^(|0) + |l))®(| 0) + |l)), (26) 

continuous measurement will eventually collapse onto one 
of the basis states and the corresponding measurement 
voltage will be distributed according to one of the four 
Gaussian distributions in Fig. 3. The smaller the relative 
phase shift of the coherent states entangled to different 
qubits levels, the longer it will take for the Gaussian dis¬ 
tributions to separate and to clearly distinguish the basis 
states. 

In contrast to the diagonal elements of the density 
matrix, which as stated above are independent of the 
frame in which the simulation is made, the coherences of 
the qubit populations will differ greatly between frames. 
Here we present an example time evolution that illus¬ 
trates the points raised in the derivation of equations of 
motion for the qubits in the lab frame in sections H G and 
HD. We assume no loss and identical cavities for simplic¬ 
ity here (we shall refer to this as the “ideal setting”), with 
simulation parameters: bi = 1: fy = ^ = 1.5MHz, = 

^ = 0.5MHz, Ai = A2 = 0 and 7^ = 72 = 0. Fig. 4 
shows the poiio off-diagonal component of the qubit den¬ 
sity matrix under a particular measurement trajectory, 
simulated using the three dynamical equations derived 
above: polaronRq. (16), compensated polaron Eqns. (18)- 
(21), and reduced Eq. (25). The top panel of the figure 
also shows the measurement pulse Adit) that produced 
the trajectory. We choose Adit) = Edit) = 0 since it is 
the appropriate drive of the two cavities in this ideal pa¬ 
rameter regime (see section HIC). As this figure shows, 
and as expected from the derivation in section HD, the 
compensated polaron and the reduced evolution equa¬ 
tions produce exactly the same results. However, both of 
these differ from the results of the polaron frame evolu- 














tion, Eq. (16), unless the photon population of both cav¬ 
ities is zero, which occurs only at t = 0 and at long times 
after the measurement pulse has ended. Another inter¬ 
esting feature is the revival of coherence produced by the 
two lab frame SMEs (polaron compensated and reduced) 
between 2 and 2.5/iS, which is after the pulse has decayed 
to zero. This is a signature of the non-Markovian nature 
of the evolution equations resulting from elimination of 
the strongly entangled cavity states. Note that a similar 
effect should be expected even in the single cavity case 
[22], and even in that simpler case, the stated lab “de¬ 
phasing rate” would be expected to be non-Markovian, 
and lead to non-decaying coherence in some parameter 
regimes. Therefore, care should be taken with the inter¬ 
pretation of parameters entering laboratory frame qubit 
SMEs derived in this fashion. 


III. CONDITIONS FOR GENERATING 
ENTANGLEMENT 


Having developed a model for the apparatus in Fig. 
1 , we now demonstrate how a sequential coherent probe 
and subsequent homodyne detection can result in a half¬ 
parity measurement of the qubits, which in turn can 
probabilistically entangle the qubits. The derivation of 
a reduced equation of motion under continuous weak 
measurement for the qubit degrees of freedom alone, 
Eq. (25), allows us to identify the exact qubit observable 
that is being monitored by the sequential probe, namely 
Re{e’‘'^{-y/Kir]ir\a + ydUi^rifc)}. In this section, we will 
specify how the parameters that can be controlled in- 
situ {e.g., frequency, amplitude and phase of drive tones 
Ad{t), Bd{t)) can be tuned such that monitoring this ob¬ 
servable can generate entanglement between the qubits 
even when the system parameters are not ideal. 

The key fact that the probabilistic entanglement 
scheme relies on is that under the half-parity measure¬ 
ment, the states jOl) and jlO) are indistinguishable. In 
that situation, starting in the initial separable state in 
Eq. (26), for entanglement to be generated by measure¬ 
ment the initial coherence between jOI) and jlO) must be 
preserved (or at least not decay substantially). For this 
to happen, the indistinguishability between these states 
must be maintained at all times, be., it is not sufficient 
that the measurement voltage be the same for both states 
merely at the final time. Specifically, we require that 
the monitored observable have the same value when the 
qubits are in state jOI) or in state jlO), or equivalently, 
for all t, = Rei?^^°^(t), which is guaranteed if 

tr [{-^yKir]ir\a{t) + joi) (oij ] 

= tr [{-y/Kiriir\a{t) + y/i^nb{t)) jlO) (lOj ] 

^ tr [{-y/Kir]ina{t) + v^nb(t))H] = 0, (27) 

where H = jOI) (Olj — jlO) (lOj. We note that derivatives 
with respect to time of the expression on the left should 


ideally also be zero (since we demanding that the condi¬ 
tion holds for all t). We will use the derivatives of Fla 
and rif, in the following and so we explicitly write the 
first derivatives of these operators here: 

na(t) = -RiHait) - iXl(Tlr\a{t) + Ad{t) 
nf,(t) = -R2nb{t) + Ki2r\ait) - ix2(Tlnb{t) 

+ Bd{t) (28) 


with Ri = Kil2 + 7i/2 -|- lA^. These equations were ob¬ 
tained by using the definitions of the operators in Eq. (11) 
and the conditional cavity state equations of motion in 
Eq. (9). The second derivatives of Fla and Flf, can be 
obtained in the same manner. 

To derive a prescription for tuning the experimental 
parameters (specifically the compensating field, Bd{t)), 
we write Fit, in terms of its first derivative using Eq. (28) 
and substitute the result into Eq. (27), to obtain 


tr 


^(Kl2Fla -|- Bd 


lib) 


K2 - iX 2 (^l 
R-^ + X2^ 


A 

) 


= 0 , 


(29) 


This is a general dynamical condition that the parame¬ 
ters in the system need to satisfy as the system evolves. 
However, this is a self-consistency equation for Bd be¬ 
cause the operator FI;, depends implicitly on Bd, and thus 
it does not provide an explicit solution for the compen¬ 
sating field. In the following, we discuss simple explicit 
solutions of Eq. (29) for three limiting, but physically 
relevant, regimes, as well as the more complex general 
solution that requires shaped pulses. 


A. Adiabatic regime 


Consider the regime where the probe fields vary very 
slowly or not at all {e.g., steady-state or continuous- 
wave measurement), specifically, this is the limit where 
Ad, Bd Ki, K 2 - In this case, we can approximate Ad{t) 
and Bd{t) as constant fields for short times and solve for 
the “adiabatic values” of the hybridized field operators 
by setting the derivatives in equations Eq. (9) to zero, 
resulting in 


nf(t) 

Uf{t) 


Ad{t){Ri - ixif^l) 

R-i^ + Xi^ 

(Ki2n°^(t) -k Bd{t)){R2 - iX2crl) 

^2^ + X2^ 


Substituting Fl^ —>• Fl^"^ in Eq. (29), dropping the lib 
term (since this is small in this regime), and then solving 
for Bd{t) in terms of the other quantities, yields 


Bf{t) = K,2Ad{t) 


(xi^2 - X 2 R 1 ) - (^2^ -b X 2 ^)Xlll ^2 
X 2 {ki^ + Xl^) 


(31) 
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FIG. 5: (Color online) Top panel: Coherence between the |01) 
and 110) states. The blue (lower), purple (middle), and red 
(upper) lines correspond to ki = {1, 5,17}MHz , respectively, 
with + 2.5 in all cases. Bottom panel: Conditional 

output held amplitudes when the qubits are in states |01) 
(Re{.B^°^^}, solid) and |10) (RejS*-^®^}, dashed), for different 
values of the cavity decay rate ki. All other parameters are 
described in the main text. 


This equation defines the value of the compen¬ 

sation held driving the second cavity, that achieves the 
desired indistinguishability of the states |01) and |10). 
Equivalently, another approach to achieving indistin¬ 
guishability, without using the compensating drive (i.e., 
setting Bd{t) =0), is to tune the frequency of the drive 
Ad- This gives a condition on Ai(<) in order to meet 
Eq. (31), and was the approach chosen in Ref. [18] due 
to its simplicity. However we note that this is only possi¬ 
ble for certain parameter ranges. In both cases, Eq. (31) 
can be met by tuning the parameter (s) in-situ to obtain 
the same measurement statistics when the qubits are in 
state 101) and in state |10). 


B. Bad cavity limit 

The adiabatic approach works best in the limit of very 
large cavity decay rates Ki, where the transient-evolution 
time periods leading up to and following the steady-state 
are short enough enough to be entirely negligible. In such 
a case, we can instead simply use a square measurement 
pulse Adi as well as for the compensation via auxiliary 
measurement drive Bd or the frequency calibration Ai. 

If we additionally assume ki,K 2 ^ Xi:X2 (i.e., take 


the bad cavity limit), Eq. (31) further simplifies to 

= - i^^Adjt) ~ ^2/2)xi + ^i^2X2/2) 

V «2 Kl^X2 

( 32 ) 

Here we have retained the time dependence for generality 
but this is typically not necessary. The simplicity of this 
approach makes it of great practical utility and as such 
it was used in Ref. [18]. Note however, that the cavity 
decay rates can only be increased up to a certain point 
imposed by physical constraints and so a small transient 
error may remain. 

Fig. 5 (upper panel) shows the degradation of coher¬ 
ence as a result of such transient populations in the cavi¬ 
ties. These transients exist because we have used Eq. (31) 
(which reduces to Eq. (32) for larger values of cavity de¬ 
cay rates) for compensation. This figure shows the decay 
of the off-diagonal element Ipoiiol while the measurement 
pulse is applied, for ^ taking on values {1,5,17}, with 
1^ = ^ -I- 2.5. ^ = 1.2MHz and all other parameters are 
the same as in Fig. 4. The total measurement time re¬ 
quired is kept approximately fixed by setting Ad = 
and employing the same pulse shape as in Fig. 4. These 
calculations employed the polaron frame dynamical equa¬ 
tion Eq. (16), with initial state iTg) and averaging over 
1500 trajectories that resulted in the desired outcome. 
The bottom panel of the figure shows the conditional out¬ 
put fields Re{i?^°^^} and Re{.B^^°^} corresponding to the 
qubits being in state |01) and |10), respectively^. Ideally, 
these conditional output fields should be identical at all 
times. However, as expected the qubit states are distin¬ 
guishable by the output fields during the pulse transients, 
since the simplified compensation prescribed by Eq. (31) 
(or Eq. (32) for the larger values of cavity decay Ki_ 2) 
is utilized. Furthermore, the greater this distinguishabil- 
ity, the greater the associated loss of coherence. Fig. 5 
shows that with larger cavity decay rates, these tran¬ 
sients become smaller; a cavity decay rate of ki = IMHz 
causes 27% loss of coherence, while in contrast a value 
Ki = 17MHz causes negligible loss. 


C. Ideal system parameter regime 

In the ideal case, when the transmission is loss¬ 
less (r]i=l) and the cavities and qubits are identical 
{ki=K 2 =k, Xi=X2 =X, Ai=A2 =0,7 i=72 =7), Eq. (31) 
reduces to 

Bdit) = -Adit). (33) 

This is equivalent to using the probe field together with 
only the reflection mode of the cavities, such that the 


® We have chosen <!> = n/2, which is shown below to be the correct 
quadrature to measure for the half parity measurement. 
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compensation fields are not needed, i.e., Ad{t)=Bd{t)=Q 
in Eq. (5). We note that the same result is obtained if 
both cavities are driven only through their input ports, 
with ^/^Bd{t) = —y^Ad{t), and e{t)=0, i.e., the probe 
field is absent. 

In this ideal parameter case it is especially simple to 
see how the sequential probe reproduces the half-parity 
measurement. Explicitly, we observe that the adiabatic 
compensation in this ideal case results in Eq. (30) reduc¬ 
ing to 

Ad{t){K/2 - ixfrl) 
k^/4 -|- 

k'^/A + 

(34) 


n 


ad,ideal 


it) = 


n 


ad,ideal 


it) = 


where we have set = 0, since the fields entering the 
input ports are zero in this case. 

The qubit observable being monitored in this ideal case 
is then 


= y^Re[e^'tAdit) 

+* X - ^ (o-z + C^z) 


ad,ideal 


)) 


Ad 


KX 

d 


i<^Wl) 


(35) 


where d = /A + x^- Choosing Adit) real and setting 
the homodyne phase 6 = '^ results in a measurement 

We note that this ideal case affords the significant 
benefit that, the transients will exactly cancel for all 
parameter values and the indistinguishability condition 
(Eq. (27)) will be met at all times. This is apparent in 
Fig. 4, where we see that in this setting, the coherence 
of the entangled state does not decrease at all during the 
periods of transient evolution of the cavities. The reason 
for this can easily be understood by looking at the equa¬ 
tions of motion for the cavity dependent qubit projectors, 
Eq. (28), in the limit of identical cavity parameters and 
using Eq. (33) for the ideal setting. Then, the condi¬ 
tion for indistinguishability for identical cavities (i.e. for 
Re{e*‘^(nf,(t) — na(t))}) differentiated twice 


tr 


(iitit) - n,(t)) 


= 0 


gives, after substituting the derivative of Eq. (28) fol¬ 
lowed by Eq. (28) itself and Eq. (33): 


tr 


x2(n,(t)-n„(t))--(n,(t)-n,(t)) H 


= 0 


From this equation it is clear that if the condition is met 
initially, it is met also at all later times. This provides 
significant motivation to make the parameters for the two 
cavities as similar as possible and to make the transmis¬ 
sion between cavities lossless, i.e., r]i = 1). 



FIG. 6: (Color online) Example of a measurement pulse Ad{t) 
and compensation pulses Edit) for two cavities with identi¬ 
cal frequencies but different transmittivities, ^=3.9MHz and 
^=3.5MHz. All other parameters are the same as in Fig. 3. 
Blue solid line: pulse Ad{t). Orange dotted line: compen¬ 
sation pulse Bd(t) when the adiabatic approximation is used 
(Eq. (31)). Red dashed line: exact solution given by Eq. (37). 


D. General dynamic condition 


Achieving the ideal setting of symmetric cavities and 
lossless transmission between them is experimentally 
challenging. In this subsection we show how to eliminate 
the detrimental effects of transients even in nonideal cases 
by shaping the measurement and compensation pulses. 
In particular, we show how the compensation held Bdit) 
can be shaped to ensure that the indistinguishability con¬ 
dition, Eq. (27), is met at all times. 

In the following, we leave Adit) unchanged (thus fully 
specifying na(t)) and shape the compensation held by 
adding a component ABdit) so that Bdit) = B^^{t) + 
ABd{t) enforces the entanglement criteria, Eq. (27), at 
all times. To derive the form of ABdit), we supplement 
Eq. (27) with its hrst and second derivative, which must 
also be equal to zero. 

The second derivative of Eq. (27) gives, upon inserting 
the derivative of Eq. (28), 


tr 



Kl2na + (k 2 -b iX2Crl)t\b 


= 0 , 


while the hrst derivative of Eq. (27) gives simply 


2K2tl' 


K 2 


= 0 . 


(36) 


Adding these to + xi) times Eq. (29) yields 


ABdit) = 

2X2 


{k2 - 1X20-1 - 


2 1 ^ 2 ^ + X2^ 


An„ 


K2 


+(i- 2 ^)n„--n, H 

K 2 K 2 


( 37 ) 
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with Ana(t) = na(t) — n“‘^(t). Eq. (37) fully speci¬ 
fies the shape parameterization of the compensation field 
Bd{t) = + A_Bd(t). Note also that the derivative 

of Ha is proportional to its deviation from the adiabatic 
value, no(t) = {ki — ixi(T],)/W\a{t). Thus, in the bad 
cavity limit (k ^ y), we find that the change in the 
compensation field relative to its adiabatic value is sim¬ 
ply proportional to the deviation of the first cavity field 
from its adiabatic value. 

This solution for the optimal shaped pulse is demon¬ 
strated for an example of two cavities with identical fre¬ 
quencies but different transmittivities in Fig. 6. Here the 
solid blue line shows the pulse the dotted orange 

line the compensation pulse Bd{t) when the adiabatic 
approximation (Eq. (31)) is used, while the dashed red 
line shows the exact solution of Eq. (37), which gives the 
optimal pulse shape that ensures the indistinguishability 
condition is met at all times. The parameters employed 
in this simulation are specified in section HE, with the 
only difference that we take unequal cavity transmittiv¬ 
ities, ^=3.9MHz and ||=3.5MHz, for this calculation. 
The exact compensation field is seen to be similar to the 
original pulse and the adiabatic approximation, with only 
slight changes during the transient periods of Ad{t). This 
ability to maintain the indistinguishability condition at 
all times by optimal shaping of the compensation pulse 
Bd{t) is very relevant for experimental situations such as 
that in Ref. [18], where the qubits could be tuned to the 
same frequencies but the cavity losses are not identical. 


IV. RESULTS FOR LOSSY TRANSMISSION 

In this section we describe results from simulating the 
dynamics of the system in Fig. 1 using the theoretical 
description developed in section II together with experi¬ 
mentally realistic parameters. We use the polaron frame 
reduced master equation, Eq. (16), with a large num¬ 
ber of random realizations of dW{t) (60000 trajectories). 
This choice of frame is made here for computational effi¬ 
ciency, since we will only study observable values at the 
end of a measurement, when the cavities are unpopu¬ 
lated and hence the polaron and lab frames coincide. For 
these simulations we use the parameters ^ = 1.2MHz, 
= l.OMHz, = 18MHz, = 16MHz, Ai = A2 = 0, 
7 i = f^,72 = and y/^Ad(t) = lOMHz. These pa¬ 
rameters are representative of the parameter regime cur¬ 
rently accessible in superconducting cavity-QED archi¬ 
tectures [18]. The compensation field, Bd{t), is chosen 
according to Eq. (31) in order to ensure the distinguisha- 
bility of the states in the single-excitation subspace, ex¬ 
cept during transients (Eq. (37) was not used in this cal¬ 
culation, since as described above, this full compensation 
requires complex pulse shaping and hence is experimen¬ 
tally more challenging). 

In a pulsed measurement setup the width of the mea¬ 
surement pulse dictates how well resolved the qubit states 
become. In Fig. 7 (top panel) we show how the state 



Normalized homodyne voltage 


FIG. 7: (Color online) Distribution of qubit state populations 
(over 60000 trajectory simulations) as a function of normal¬ 
ized homodyne voltage and measurement pulse width. The 
initial state is given in Eq. (26). Black (low voltage) repre¬ 
sents 100) population, blue (high voltage) jll), red (upper, 
medium voltage) :^(|01) -1- jlO)), and green (lower, medium 
voltage) represents ■^(jOl) — jlO)). The top panel plots rel¬ 
ative frequencies of each of the state populations, while the 
bottom panel plots normalized populations conditioned on the 
measurement voltage value from the x-axis. The bottom plot 
is essentially a slice of the top panel taken at a measurement 
pulse duration of 1/rs. Simulation parameters are specified in 
section IV. 


populations are distributed as a function of the (nor¬ 
malized) homodyne voltage and the measurement pulse 
width, over the 60000 simulated trajectories. The bot¬ 
tom panel shows a slice at measurement pulse width of 
I/iS, where the populations are normalized (such that 
the total population over these four orthogonal state is 
fixed to be one) at each homodyne voltage value. We see 
that for short pulse widths little information is carried 
out of the cavities and all qubit states are equally likely 
(low SNR). As a result of employing the compensation 
pulse Bd{t) we see that, for all pulse widths, the |0I) and 
110 ) states are indistinguishable by the homodyne voltage 
value. In addition, for pulse widths greater than ^ I/rs, 
we find that the homodyne voltages concentrate around 
the center value, predicting that primarily the single ex¬ 
citation subspace is populated under these conditions. 
The presence of the ^ |0I) -I- |I0) state indicates that the 
indistinguishability condition has been satisfied between 
|0I) and |I0). However, even if this is the case, other 
sources of dephasing could result in a mixed state with 
no entanglement for which the homodyne voltage would 
also be concentrated around the center value. For this 
reason we also plot the population of the antisymmet- 
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ric state ^(|01) — |10)) (green line), presence of which 
would indicate a mixture with reduced or no entangle¬ 
ment, depending on the relative value of this state. We 
see that the antisymmetric state does not contribute for 
measurement pulse widths ~ 1/rs, but that its popula¬ 
tion increases as the measurement pulse width increases, 
as a result of intrinsic dephasing becoming significant at 
longer times. Fig. 7 therefore shows that there is a trade¬ 
off between achieving high SNR with long measurement 
pulses and compensation for indistinguishability, and re¬ 
stricting the pulse duration to avoid intrinsic dephasing 
at longer timescales. 

To explore the influence of the two primary detrimen¬ 
tal effects to achieving entanglement between qubits in 
separate cavities, namely loss between cavities (1 — rji) 
and measurement inefficiency (1 — 77^), we quantify the 
maximum achievable entanglement, quantified here by 
the concurrence [27], for a range of these parameters in 
Fig. 8. The maximal achievable concurrence (maximized 
over the measurement pulse width, over a range 0.1—4/rs, 
and averaged over all trajectories) is plotted as a function 
of both loss of photons between the cavities {rji in dB) 
and measurement efficiency rym- In these simulations the 
indistinguishability criteria is enforced by the adiabatic 
compensation pulse, Eq. (31). The transmission losses 
between the two cavities lead to dephasing of the first 
qubit, and hence degrade the entangled state formed with 
the second qubit. On the other hand, decreased measure¬ 
ment efficiency does not by itself lead to lower coherence, 
but instead necessitates longer measurement pulse widths 
or larger probe field amplitudes (ie., we need to use more 
photons to obtain the same amount of information). Nev¬ 
ertheless, low measurement efficiency combined with loss 
between cavities is detrimental, because even though one 
can increase the number of photons used to probe the 
qubits, this results in more photons being lost between 
cavities and thus in greater dephasing. Therefore the 
combination of low measurement efficiency together with 
large transmission loss between cavities is the most un¬ 
favorable situation for generating qubit entanglement. 


V. POPULATION TRANSFER IN THE 
POLARON FRAME 

The reduced model developed above and the analy¬ 
sis that followed are only valid when the original master 
equation, Eq. (3), contains no terms that do not commute 
with the polaron transformation. Two such terms that 
one might like to include in an extended model are qubit 
driving and population relaxation (T1 process). In this 
section we show how to perturbatively incorporate these 
effects into the qubit-only reduced master equations in 
Eq. (16) or Eq. (25). 



Measurement efficiency 


FIG. 8: (Color online) Maximum concurrence versus measure¬ 
ment efficiency (?7m) and photon loss in the channel between 
the cavities {1 — rji, in dB). Simulation parameters are spec¬ 
ified in section IV. Entanglement maximum is calculated by 
propagating the simulation for different measurement pulse 
widths and selecting the maximal concurrence achieved for 
the given parameters averaged over all trajectories. 

A. Qubit driving 

Resonant drives on the qubits are described by addi¬ 
tion of the Hamiltonian term 

Hd = ^i{t)(Tl+n2(t)(Tl. (38) 

We will study instances where the magnitude of these 
new driving terms is small, and therefore they can be 
treated perturbatively. 

As a result of the hybridization of the qubits and cav¬ 
ities in the polaron frame, qubit transitions will result in 
the cavity also being driven. This can be seen from the 
effect of the two-cavity polaron transformation, Eq. (10) 
on the cr!_ operators: 

{crlf = D2[nb]Di[na](TlDl[na]Dl[nk] 

(^2_)P = D2[ni,]<TlDl[n,]. (39) 

Eor simplicity, we consider an adiabatic probe pulse and 
compensation pulse. Then using the parameterization of 
Eq. (30) we obtain to first order in the dressed 

cr!_ operators 

{cr\)^ « cr\{l — Adi'a + C,2^)iii + A*d{a^-\- 

« <T^(l-(Rd + A,Ci)Af2b + (H2+A5Cr)Af;bI), 

(40) 

where /i* = 2xj/((Ki/2-|-zAj)2 -bxf) and Ci = Ki2(Ki/2-l- 
iA,)/((Kj/2-biAi)2 -bxD- 

These dressed operators contain both qubit and cavity 
operators and thus in the polaron frame, qubit driving 
will also lead to cavity driving and damping. Since the 
cavities are detuned from the qubits the counter-rotating 
terms in the above expansion can be dropped in the 
RWA, and one obtains cavity sideband transitions involv¬ 
ing the qubits (notably, driving the first qubit can result 





13 


in a sideband transition in the second cavity due to the 
interconnection). Further, in the strong cavity damping 
limit, when ^ > the qubit driving Hamiltonian 

in the polaron frame, can be approximated as 




Hicr^ - 


+ Hixicri) 


( ki /2 + iAi)2 - yi2 _ 

+ HiXicrj,) 

{K 2 / 2 +lA 2 )^-Xl^-m^ 

2 ^ 2 ^ al + ^ 2 X 20-1 


+ H 2 <t — 2 A'^ 


(K 2/2 + zA2)2 - X2^ - 4H2' 


r(41) 


with X=Bd^i 2 + Ad^i 2 C,i- The cavity coupling induces 
a drive-dependent energy shift of the qubit, or alterna¬ 
tively, tilts the drive axis. This approximate drive Hamil¬ 
tonian can be used with the qubit-only master equation, 
Eq. (16). However, effects outside of the qubit subspace 
will also remain in the polaron frame, such as sideband 
heating of the cavity. It is evident from Eq. (41) that the 
same control pulses will result in different angles of rota¬ 
tion, when viewed in the laboratory and polaron frames, 
and will only agree in the limit of zero photons in the cav¬ 
ities. In the dispersive limit and in steady state, the axes 
of rotation Eq. (40) are constant and it is then straight¬ 
forward to specify the X and Y rotation angles to achieve 
the desired rotation in either frame. 

In addition to this tilting of the qubit drive axis, the 
simultaneous application of a measurement tone and a 
qubit drive can lead to measurement-induced suppression 
of coherent oscillation (the Zeno effect), as noted for the 
single qubit case in Ref. [22]. 


equation, Eq. (3), by the addition of the Lindblad terms 

2 

CrQ = Y,xlV[cTl]g, (42) 

i=l 

where 7 ® are the intrinsic relaxation rates. The ap¬ 
proximation Eq. (40) to the polaron frame form for (T*_ 
Eq. (39) illustrates the problems that will arise when at¬ 
tempting to carry out the derivation of reduced master 
equations as in section H in the presence of this relaxation 
term. Specihcally, the presence of field excitation terms 
in the polaron frame mean that the cavity states are now 
no longer in the vacuum at all times in this frame. Phys¬ 
ically, this reflects the fact that the relaxation process 
can create cavity photons which are not compensated 
by the polaron frame transformation. This means that 
Eq. (16) is no longer a valid equation of motion for the 
qubit degrees of freedom in the polaron frame. Further¬ 
more, Eq. (17) is not a valid ansatz for the state of the 
qubit and cavity degrees of freedom. However, as we will 
show below, when the relaxation is slow and in the bad 
cavity limit, one can return to the lab frame and derive 
a valid reduced equation of motion in this frame that 
incorporates the qubit relaxation. 

We begin with the more general representation of the 
state of the system in the polaron frame, that takes into 
account that the cavity states are no longer always the 
vaccum in this frame: 

rn^nbm^mb,ijkl{t) \ij) {kl\®\nanb) (TOaWbl , 

rLaTibrnamb 

ijkl 


B. Qubit relaxation 

For long measurement pulses or continuous wave mea¬ 
surement, incorporating the effects of qubit relaxation 
(Ti process) can become important. Similar to intrinsic 
dephasing, relaxation is incorporated into the full master 


where the indices i,j, k, I run over 0 and 1 , while the other 
indices (that index photon number states) run from 0 to 
00 . Carrying out the same reduction to a representation 
of the state of the qubits in the lab frame made in section 
HC (Eq. (19)), one finds that the mapping between the 
diagonal and off-diagonal terms in the polaron frame and 
the reduced lab frame is: 


Pijij{t) — 

Tlaflb 

(43) 

Pijkl{p^ — 

rn,ntm^mb,ijki{t) {mamb,kl\Dl[r\a{t)]Dl[r\b{t)]Di[na{t)]D2[r\b{t)] \nanb,ij) 



YlaflbrnaTTlb 


= 

^ ^ ^ria.'ribTnarnb^ijkli^^ ■ 

(44) 


naTihmarrib 


In the presence of qubit relaxation, calculating the time 
derivative of these off-diagonal elements is more involved 
since we cannot assume that roooo,ijki is the only nonzero 
value in this expansion. 


Let us focus on the effect of the relaxation term 
Eq. (42) alone. Incorporation of this term adds the fol¬ 
lowing to the derivatives of the diagonal terms (ignoring 
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all the other terms in the dynamical equation): 


710 ^ 5 , 0000(0 


Tr ^710 715770,715,1010 (^) 
+ 7 r^ 71 o 7 l 571 on 5 , 010 l(^) 


^7lo7l5 71o7l5 ,0101 (0 

^TIq 77-5 Tlo 715,1010(0 

^na7l571o7l5,llll (0 


7r ^Tlo 715 Tlo 715,1111 (^) 

7r ^Tio 715710 715,0101 (^) 

7r ^7lo7l571a7l5,101o(^) 

+7r^7i 0715710 715,1111 (^) 

~ilr “1“ 7r)^n.on6"ani,,llll(0- 


the lab frame reduced SME can be approximated by the 
addition of the following Lindblad term to Eq. (25) 

Crp{t) = 7,(X>[cr!_]p(t), (47) 

i=l,2 

and hence qubit relaxation carries through unaltered to 
the reduced SME in the lab dream. 


VI. SUMMARY AND CONCLUSIONS 


When the sum over Ua and ni, prescribed in Eq. (43) is 
performed, one gets equations of motion for Pijij (t) that 
are consistent with a qubit relaxation process in the re¬ 
duced lab frame. However, the effect on the off-diagonal 
components is not as straightforward. To illustrate this, 
we shall focus on a single off-diagonal element and cal¬ 
culate the contribution of the relaxation term ^^'D[cr^]g 
to the evolution of poooi(0- Explicitly, we obtain (again 
ignoring all the other terms in the dynamical equation), 


POOOl(0 — ^ '^n„ra6m„mt,OOOl(0 

7la7l57na77i5 

“ 7r ^ ^ '^7la7757na7n5,1011 (^) 

TianhmaTnh 

7r ^ ^ ^71(1,7757710,7715,1011^71a ,771^ 

TlaTLimaTni, 

Jr / j ' 71a7l571(i7715,1011'^ 

Tlanbmf, 


The summation in this expression is difficult to perform 
exactly. However, we may simplify this by assuming the 
bad cavity limit and using the adiabatic values for the 
cavity fields given in Eq. (30). In this situation. 




-»2x2-Bd 

fil + xl 


< 1 , 


(45) 


where the inequality in a consequence of the bad cavity 
limit. Similarly, in this limit we have « 

1. Now the matrix element (TOf,| i72[A] jn^) is propor¬ 
tional to the Laguerre polynomial L”|j“™'>(A), and is 
peaked around na = ni, and zero everywhere else for 
small X [28]. Therefore we approximate (mbj — 

\n}y) « and in this bad cavity limit (again 

for adiabatic values of the cavity fields), we find 


POOOl(i) ~ 7r f’ria.mnanb,— 7rPl01l(i)- (46) 

71a7l5 


Calculating the equations of motion for the other off- 
diagonal elements under the same approximations, we 
find that the total contribution of qubit relaxation to 


We have derived a framework for describing joint dis¬ 
persive measurement of qubits in separate cavities and 
shown how these measurements may be used to engi¬ 
neer entanglement between pairs of such qubits. The de¬ 
scription shows how the populations of the qubit levels 
evolve diffusively as a function of a cascaded measure¬ 
ment and how calibrating the amplitude or frequency of 
compensation field(s) enables preservation of the coher¬ 
ence between the single-excitation states during the mea¬ 
surement. We derived two alternative stochastic master 
equation approaches for the dynamical description of the 
qubit density matrix, one based on a polaron representa¬ 
tion that describes dynamics in a dressed frame and the 
other a reduced master equation that described dynam¬ 
ics in the bare laboratory frame. We derived static and 
dynamic entanglement conditions in the energy basis of 
the system, i.e., conditions on the physical parameters 
that guarantee the monitoring of the output field will 
ensure entanglement of the two qubits, and showed that 
these conditions give a simple prescription for ensuring 
this indistinguishability at all times. We discussed simple 
solutions of these conditions in three physically relevant 
limiting regimes of adiabatic probes, bad cavities, and 
an ideal setting of identical cavities with lossless trans¬ 
mission. We further showed how to achieve the indis¬ 
tinguishability condition at all times by optimal shaping 
of the compensation field Bd{t). The theoretical analy¬ 
sis was applied to realistic experimental situations with 
extensive simulation of the different stochastic master 
equations in polaron and in laboratory frames. The sim¬ 
ulations show that the entanglement achieved with this 
procedure is tolerant to significant imperfections in the 
measurement efficiency and, to a lesser extent, to the 
presence of loss along the probe field path, e.g., arising 
from circulators to ensure unidirectionality. A detailed 
comparison was made between the simulations of the dif¬ 
ferent qubit stochastic master equations in the polaron 
and laboratory frame. A single trajectory analysis and 
related simulations revealed that the intra-cavity fields 
provide a non-Markovian environment resulting in sup¬ 
pression and revival of coherence of the bare laboratory 
qubit states, indicating that use of Markovian models for 
reduced qubit equations of motion in the lab frame is not 
always accurate. 

The theoretical formulation presented here provides a 
first-principles description of the remote probabilistic en¬ 
tanglement achieved experimentally in Ref. [18]. In addi- 





15 


tion, these results also motivate other schemes to achieve 
entanglement of superconducting qubits. Of particular 
interest is the extension to continuous wave measure¬ 
ments rather than the pulsed measurements that were de¬ 
scribed here. Continuous wave measurements are needed 
for longer and more complicated applications such as er¬ 
ror correction protocols and feedback control in order to 
stabilize the entangled state in the presence of dephasing 
and relaxation processes. We have shown in this work 
that in the weak driving limit the main effect of simulta¬ 
neous measurement and coherent control is a tilting of the 
axis of rotation that can be compensated (and a suppres¬ 
sion of the rotation due to the measurement in the strong 
measurement limit, as was shown in Ref. [22]). We have 
also shown that in the bad cavity limit phenomenologi¬ 
cal qubit relaxation terms simply carry over to the lab 
frame reduced master equation for the qubit degrees of 
freedom. This motivates introducing measurement-based 
feedback to create multi-qubit entanglement determinis¬ 
tically, which will be the topic of a future publication. 
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Appendix A: The SLH representation of quantnm 
networks 

Building off work by Gardiner [19] and Carmichael [20] 
on cascaded quantum optical systems, Gough and James 
have constructed a general formalism for modeling net¬ 
works of quantum systems connected by bosonic fields. 
The utility of the approach is that it enables description 
of the dynamics of complex networks of modular compo¬ 
nents using simple composition rules. This formalism is 
primarily developed in Refs. [23, 24] and we summarize 
the main results in this Appendix. 

The basis of the SLH modeling approach is to decom¬ 
pose a network into localized components with arbitrary 
degrees of freedom that are connected via freely propa¬ 
gating unidirectional broadband fields. This allows one 
to eliminate the fields propagating between components 


to arrive at an effective description of the system just in 
terms of the localized degrees of freedom and how they 
are connected together. 

The starting point for this formalism is the Hudson- 
Parthasarathy quantum stochastic differential equation 
(QSDE) accounting for time evolution of the unitary op¬ 
erator, U{t), describing coupled evolution of the system 
and field degrees of freedom [29]: 

dU{t) = {{S-l)dA{t)+LdB'<{t)-L'<SdB{t) 

-{^L'<L + iH)dt\uit), (Al) 


where B{t) and are integrated versions of the freely 
propagating bosonic fields linearly interacting with the 
system at an interface or “port” (these could be output 
fields from another system): 


B(t) = f b(s)ds, B\t) = f {s)ds, (A2) 
Jo Jo 


with [6(t),&l(s)] = S{t — s). This commutation relation 
defines the bosonic fields as rather singular objects, and 
hence the increments, dB{t) = B{t + dt) — B{t) (and sim¬ 
ilarly of dB\t)), are operator valued stochastic variables 
that are analogous to Ito increments. Finally, A{t) is a 
quantum stochastic process that corresponds to the ob¬ 
servable counting the number of quanta in the bosonic 
field that have interacted with the system up to time t: 


m 



b\s)b{s)ds 


(A3) 


The other components of Eq. (Al), the system opera¬ 
tors S', L, and H, describe the nature of the interaction 
between the system and propagating field at the inter¬ 
face. S describes the impact on system when photons are 
scattered between ports (this component is most interest¬ 
ing when we consider systems with multiple ports, as we 
shall below), L is the system operator that is directly 
and linearly coupled to the field, and H is the system 
Hamiltonian that accounts for dynamics that does not 
involve interaction with the field b{t). These components 
are often grouped together into a triple G = {S,L,H), 
which is sufficient to completely characterize the system 
evolution. 

The generalization of Eq. (Al) to the case where the 
system has multiple ports, with independent fields at 
each port interacting with system, is: 


dU{t) = - Sjk)dAjk{t) 


jk 


SjkdBkit) 


jk 


+^H)dt [ U{t), 


(A4) 
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where Sjk describes the effect on the system of a photon 
scattering from port j to k, and Lj is the system operator 
coupled to the field at port j. In this multi-port case, we 
still describe this localized system with an SLH triple, 
but with the S and L now being matrix or vector valued: 



Note that the components of the L vector are themselves 
operators. 

The key advantage of the SLH formalism is that one 
can easily construct effective descriptions of arbitrar¬ 
ily connected networks of localized components, each of 
which is represented by a triple: G = {S, L, H). Connect¬ 
ing two components in series, parallel, or in feedback re¬ 
sults in another system represented by another SLH triple 
whose matrices can be derived by simple algebraic rules 
[23] . For example, consider connecting two localized sys¬ 
tems in series, where the outputs from Gi = {Si, Li, Hi) 
and connected to the inputs of G 2 = {S 2 , L 2 ,H 2 ), where 
for simplicity we assume that the number of input ports 
that G 2 has is the same as the number of output ports 


that Gi has. The resulting system is represented as 
G 3 = G 2 <\Gi 

= {Sz = S2Si,L^ = S2Li-\-L2, 

iLs = iJi-biLa+Im|45'2Li}). (A6) 

Similarly, if one connects Gi and G 2 in parallel (concate¬ 
nates them), the resulting system is represented as 

G3 = GsfflGi 



H:i = Hi+H2). (A7) 

See Refs. [21, 23] for more details on more complex com¬ 
position rules. 

Appendix B: SLH triple for cascaded cavity 
apparatus 

The SLH triple that results from performing the con¬ 
catenation and series products in Eq. (2) is 


r 

-iy/l - Pi 

0 

0 1 


' VI - Pi{-e{t) + y/Kia) ■ 


^/m 

0 

0 


y/rfie{t) - y/Kipia + 

0 

0 

-1 

0 

1 

— Ad{t) -b 

L 0 

0 

0 

-ij 


L -Bd{t) -b J 


Aia'l'a -b xia'l'acr^ -b A2b^b -b - aA2(t)) -b ^(b'l'Hd(t) 


bB:(t))-^(atb-bta))j(Bl) 


For any system represented by an SLH triple {S,L,H), 
the corresponding master equation that describes the dy¬ 
namics of the internal states in the model, represented by 
the density matrix p, is: 

= -i[H, p]+'^LupLl - ^LlLkp-^pLlLk, (B2) 

k 

where Lk are the (operator-valued) elements of L. No¬ 
tice that the scattering matrix S does not influence the 
internal dynamics of the system, it only determines the 


relation between the input and output fields. 

Finally, using the fact that the evolution of a density 
matrix under Eq. (B2) is invariant under the transforma¬ 
tions 

Lik Lk + O! 

H ^ H - "-{Lka* - Lla), 

for a G C, we can remove the terms proportional to iden¬ 
tity in the L vector to obtain the equivalent SLH triple 


^/m -VI - 0 0 ■ 


yVa) 

-VI - m vV 0 0 


-y/Kipia + V^b 

0 0-10 



L 0 0 0 -ij 


L V^b J 


Aia'l'a + xia'I'acr^-b A 2 b^b-b X 2 b'^b(T^-b i(a'l'Ad(t) - aH2(t)) + i(b^Hd(t) - bH2(t)) - ^^(a'^b - b^a)) j (B3) 
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After addition to phenomenological dephasing terms with 
rate parameters = 1,2, this SLH triple corresponds 
to the full master equation given in Eq. (3). Each element 


of L corresponds to an “output” port that leaks photons, 
and the second port is the only one that is monitored 
(see Figs. 1 and 2). 
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